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We extend the formalisms developed in Gair et al. [T] and Cornish and van Haasteren to create 
maps of gravitational-wave backgrounds using a network of ground-based laser interferometers. We 
show that in contrast to pulsar timing arrays, which are insensitive to half of the gravitational-wave 
sky (the curl modes), a network of ground-based interferometers is sensitive to both the gradient 
and curl components of the background. The spatial separation of a network of interferometers, 
or of a single interferometer at different times during its rotational and orbital motion around the 
Sun, allows for recovery of both components. We derive expressions for the response functions of a 
laser interferometer in the small-antenna limit, and use these expressions to calculate the overlap 
reduction function for a pair of interferometers. We also construct maximum-likelihood estimates of 
the -I- and x-polarization modes of the gravitational-wave sky in terms of the response matrix for a 
network of ground-based interferometers, evaluated at discrete times during Earth’s rotational and 
orbital motion around the Sun. We demonstrate the feasibility of this approach for some simple 
simulated backgrounds (a single point source and spatially-extended distributions having only grad 
or curl components), calculating maximum-likelihood sky maps and uncertainty maps based on 
the (pseudo)inverse of the response matrix. The distinction between this approach and standard 
methods for mapping gravitational-wave power is also discussed. 

PACS numbers: 04.80.Nn, 04.30.Db, 07.05.Kf, 95.55.Ym 


I. INTRODUCTION 

Searches for anisotropic gravitational-wave back¬ 
grounds have typically been formulated in terms of 
the distribution of gravitational-wave power on the sky 
(see, e.g., [SHSl)- The basic idea underlying these ap¬ 
proaches is to use to cross-correlation measurements from 
two or more detectors to estimate the power in the 
gravitational-wave background as a function of sky po¬ 
sition. For a network of ground-based laser interferome¬ 
ters like LIGO [TU], Virgo [TT], etc., or space-based inter¬ 
ferometers like LISA [T2], eLISA [T31 [II] or BBO [T5] . 
the motion of the detectors modulates the correlated 
gravitational-wave signal at harmonics of the Earth’s 
daily rotational motion, or the spacecrafts’ yearly orbital 


motion around the Sun. The time-varying signal carries 
information about the multipole moments characterizing 
the anisotropic distribution of power, which can be esti¬ 
mated using, e.g., maximum-likelihood methods [5]. 

Recent papers by Gair et al. [T] and Gornish and 
van Haasteren [5] describe an alternative approach 
for mapping the gravitational-wave sky, which can be 
used to recover both the amplitude and phase of the 
gravitational-wave signal at each sky position. The anal¬ 
ysis in P] is cast in terms of the traditional plus and cross 
polarization components {h+{f,k), hxif,k)}, while in 
[T] the metric perturbations are decomposed in terms of 
spin-weighted or tensor (gradient and curl) spherical har¬ 
monics ’ ^0m)ab (fc)}. This latter decomposi- 
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tion is similar to that used to characterize the polariza¬ 
tion of the cosmic microwave background m. taking into 
account the symmetric, transverse-traceless nature of the 
metric perturbations habif, k). 

Although the formalisms introduced in D m are 
general, they were applied specifically to the case of 
gravitational-wave searches using pulsar timing arrays. 
In contrast to the case of ground-based interferometers 
on a rotating, orbiting Earth or the orbiting LISA/eLISA 
spacecraft, a pulsar timing array operates effectively as 
a static galactic-scale gravitational-wave detector 
with each Earth-pulsar line-of-sight being the equiva¬ 
lent of a one-way, one-arm interferometer with a com¬ 
mon endpoint at the solar system barycentre (SSB). This 
is because the frequency range for pulsar timing mea¬ 
surements is such that the displacement of a detector on 
Earth from the SSB is much smaller than the wavelength 
of the relevant gravitational waves, and hence the detec¬ 
tor effectively resides at the SSB. In addition, the times 
of arrival of the pulses from a given set of pulsars at a 
radio telescope are typically converted to arrival times at 
the SSB, which is a fixed origin. A consequence of the 
static nature of pulsar-timing measurements is that a pul¬ 
sar timing array will never be able to observe half of the 
gravitational-wave sky regardless of how many pulsars 
are included in the array [1] . In particular, the response 
of a pulsar timing array to curl modes is identically zero, 
as the gravitational contribution of such modes to the 
timing residual equals zero when integrated over the sky. 

In this paper, we extend the formalisms developed in 
mm to the case of ground-based interferometers like 
LIGO, Virgo, etc. For simplicity we work in the small- 
antenna (or long-wavelength) limit, which is appropriate 
for such detectors. We show that in contrast to pulsar 
timing arrays, a network of ground-based interferometers 
is sensitive to both the gradient and curl modes of the 
background. The fact that the spatial separation of the 
detectors is the same order or greater than the radiation 
wavelength is sufficient to allow for recovery of both types 
of mode. We demonstrate this by: (i) explicitly deriving 
analytic expressions for the gradient and curl response 
functions of a laser interferometer, and (ii) constructing 
maximum-likelihood estimates of the gravitational-wave 
sky for different types of simulated backgrounds. The re¬ 
construction of the sky maps is based on singular value 
decomposition (SVD) of the whitened response matrix 
R = USVf, which maps the modes of the gravitational- 
wave background to the response of the individual inter¬ 
ferometers, evaluated at discrete times during Earth’s ro¬ 
tational and orbital motion around the sun. The columns 
of U and V corresponding to the non-zero singular val¬ 
ues of S have the interpretation of response range vec¬ 
tors and sky map basis vectors, in terms of which the 
measured response and the maximum-likelihood sky map 
can be written as linear combinations [2]. The recovered 
sky maps can be calculated in terms of either a pixel- 
based parametrisation, {h+{f,kn),hxif,kn)}, where n 
labels the pixels on the sky, or in terms of the gradi¬ 
ent and curl spherical harmonic components 


®pm) (/)}’ where (Im) labels the various multipole modes, 
up to some maximum value ^max- 

The organization of the rest of the paper is as follows: 
In Sec. [nj we summarize key formulae related to the 
tensor spherical harmonic decomposition approach de¬ 
scribed in [1]. We derive, in Sec. Ill analytic expressions 
for the gradient and curl response functions of a laser in¬ 
terferometer in the small-antenna limit. (Details of the 
derivation are given in Appendix]^) In Sec. |IV[ we show 
that we can recover the standard overlap reduction func¬ 
tions using the analytic expressions for the response func¬ 
tions derived in Sec. |III[ and in Sec. |V] we compare the 
effects of Earth’s rotational and orbital motion on the sky 
reconstruction. We describe the map-making formalism 
in Sec. |VI| and demonstrate that a network of ground- 
based inteferometers can recover both the grad and curl 
components of a gravitational-wave background, by con¬ 
structing maximum-likelihood sky maps for some simple 
simulations. We conclude in Sec. IVIII with a brief sum¬ 
mary and discussion of the results, listing a few mod¬ 
ifications that might be needed in order to apply this 
formalism to real data. 


II. BRIEF REVIEW OF TENSOR SPHERICAL 
HARMONIC DECOMPOSITION 


For completeness, we summarize in this section some 
key formulae from the tensor spherical harmonic decom¬ 
position method described in Interested readers are 
referred to that paper for more details. 

Any transverse-traceless tensor field on the sky can be 
decomposed as a superposition of gradients and curls of 
spherical harmonics: 


^{lm)abik) Ni Y(^ixa)-ab{k) 2 

^{lm)abik) 2 ^(lm)-,ac{k')e b + 


( 1 ) 


where semi-colon denotes covariant derivative, gab is the 
metric tensor on the sphere, Cab is the Levi-Civita anti¬ 
symmetric tensor 


Cab — 



and Ni is a normalization constant 


Ni 


2(1-2)! 

JTW' 


( 2 ) 


(3) 


In terms of the standard polarization tensors on the sky 

eab(^) =&aOb- 

^abO^) =^a^b + iJb, 


we have 
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yifm)atCk) = ^ (fc) 

y^m)abik) = 


^{i7n){k)e^^ik) - X(,^)(fc)e+(fc) 
where W(^im)ik) and X(^i^'^{k) can be written in terms of associated Legendre polynomials 


( 5 ) 


as 




+ m)^P,™l(cos0)U*™^ (6) 

Sin 0 ' 




+ [(;-1) cos 0Pr(cos0)-(I+ m)P,™i(cos 0)] I e™^ 


and are related to spin-2 spherical harmonics [iHiiin] through the equation 


TV; r - . 

±2^(lm){k) ^y(lm){k') iX(jjj^'^{k^ 


(7) 


( 8 ) 


In terms of the grad and curl spherical harmonics, a general gravitational-wave background can be decomposed as 

hab{t,x)= r df [ d2L!^h,,(/,fc)p2./(t-fe..7c)^ (9) 

J — oo J S'^ 


where 


OO I 

habif, k) = '^Yl [«°m)(/)^(L)a&(^) + aflni)i.f)Y^m)ab(k) 

1—2 m——l 


( 10 ) 


The mode functions related to 

the more traditional -I- and x polarisation components 
h+{f, k), hx if, k) defined by 


Kbif, k) = h+if, fc)e+(fc) + hx if, k)e^,Ck), (11) 


III. RESPONSE FUNCTION CALCULATIONS 


The frequency-domain response of a laser interferome¬ 
ter to a gravitational-wave background is given by 

f(/)= f d20^^P^(/,fc)/i^(/,ft), (14) 

Js^ ^ 


Via 


h+if,k) afi^)ifW(in,)ik) - afi^)if)X^i^~^ik) 


(Im) 


2 1 


hxif,k) ^fim) (/)^(/m) ik) + afi^) (/)kL(;™) (ft) 


(Im) 


2 1 


and 


( 12 ) 


4nrif)=Nl[h+if, ft)lU()^)(ft)+ftx(/, ft)^(*;™)(ft) 

4m)if)=NlJ^f^k [kxif, ft)l^f;4fc)-M/, k)Xllm)ik) 

where we use the shorthand J2{im) — X)m=-r 

Note that all the summations over I start at ^ = 2. 


where A = {-b, x} and 

k) = ^effc(fc)(u“u^ - y<^yb)e-^-2^fk-So/c (^5) 

Here u, v are unit vectors along the arms of the inter¬ 
ferometer and xq is the position vector of the vertex of 
the interferometer at the time t when the measurement 
is made. The above expression for R^if,k) is valid in 
the small-antenna limit, which is appropriate for ground- 
based interferometers like LIGO, Virgo, etc. (see e.g., 
[21)) for this discussion in the pulsar timing context). The 
length of the interferometer arms do not enter the expres¬ 
sions for the response functions in this limit. 

Alternatively, if we expand the metric pertur¬ 
bations in terms of the gradient and curl modes 
{“(/m)(/)’ ®(i'm)(/)}’ response can be written as 

^(/) = H 

{Im) P 
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where P = {G, C} and 


R' 




R 


(Im)if) = Y 



R+if, k)W^in.){k) + R^if, k)Xi^im){k) 

[ d'llfc 

Is- 

R^ if, k)W^im)Ck) - R^if, fc)^(M W 


(17) 


These integrals were evaluated in Appendix D of [T] for the case of a reference frame which has Xq = 0: 

47r /T 


R 

R^ 




(18) 




Note that the curl response is identically zero in this frame, while the gradient response is independent of frequency 
and is non-zero only for the quadrupole components, 1 = 2. These results are qualitatively similar to those for a 
pulsar timing array, which also have zero curl response, and a frequency-independent grad response proportional to 
Yim(u), where u points in the direction to the pulsar. In what follows, we will use the notation 


47r /1 

Fm{u,v) = —'J - \Y 2 m{u) - l 2 m(D)] 


(19) 


to denote the particular combination of spherical harmonics that appear in Eq. (18). Since only the quadrupole 
response is non-zero, the index m on Fm is restricted to have values m = 0, ±1, ±2. 

For a single static interferometer, there is no loss in generality in choosing a reference frame with the vertex of the 
interferometer located at the origin, with the response functions given as above. But for a network of interferometers 
attached to a rotating and orbiting Earth, such a coordinate choice is no longer possible. If an interferometer is 
displaced from the origin by a?o, one can show that 

2 1+2 L 

R{im)if)= Rm'{u,v)‘iTr{-i)^jL{a)YlM{xo) 

m'^-2 L^l-2 M^-L 
(- 1 )” 


{2-2 + l){2l + l){2L + l) f 2 I L 
47r V ^ ^ 


2 I L 
2-2 0 


1+2 


m'— — 2 L—l — 2 M— — L 

(-1)™' l{2-2 + l){2l + l){2L + l) 


2i 


47r 


I 


L 

M 


2 I L 
2-2 0 


^ [(-l)'+^-bl] 


[(_l)*+i _ 1] 


( 20 ) 


where a = 2TTf\xo\/c and jh^cx) are spherical Bessel func¬ 
tions of order L. These expressions are derived in Ap¬ 
pendix]^ The two expressions in parentheses ( ) for each 
response function are Wigner 3-j symbols (see for exam¬ 
ple m, 122]). Note that, in general, the curl response is 
now non-zero, in contrast to the static single interferom¬ 
eter case. In addition, both response functions depend 
on frequency via the quantity a, which has the physical 
interpretation of being 27r times the number of radia¬ 
tion wavelengths between the origin and the vertex of 
the interferometer. Since the coordinate system for the 


response functions in (18) was not chosen in any par¬ 


ticular orientation relative to the unit vectors u and v, 


Eqs. (20) are valid in an arbitrary translated and rotated 
coordinate system, provided we use the angles for u, v, 
and xo as calculated in the rotated frame. 


IV. RECOVERY OF STANDARD OVERLAP 
REDUCTION FUNCTIONS 


Given the above expressions for 7?p^)(/) one can cal¬ 
culate the overlap reduction function ri 2 (/) for a pair of 
interferometers for gravitational-wave backgrounds hav¬ 
ing different statistical properties. For example, for a 
statistically isotropic background with Gp*^ = Gp'" = Ci 
and Gp*^ = 0 = Gp*^, it was shown in [T] that 


Tl2{f) ^ 12,lif)'. 


( 21 ) 


Z=2 
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where 


i 

ri2X/)= E RHim)if)R2iU(fy (22) 

m——l 

This is most-easily evaluated in a reference frame in 
which the vertex of one interferometer is located at 
the origin, and the vertex of the second interferome¬ 
ter is located along the z-axis of this frame. This so- 
called computational frame is related to the cosmic ref¬ 
erence frame located at the solar system barycentre via 
a translation by the position vector Xi and a rotation 
TZ = Tiy{Oo)TZzi4>o) such that Ax = X 2 — xi is directed 
along the z-axis. (Here {6o,(j)o) are the polar and az¬ 
imuthal angles of Ax with respect to the cosmic frame.) 
In the computational frame 

Riiim)if) = S^^5i2F^{TZui,TZvi), (23) 


where the arguments of are the polar and azimuthal 
angles of ui and vi with respect to the computational 
frame. This form for .Rf(;„)(/) implies 

2 

ri2.i(/) = fc E ^’n^(7^Ul,7^^;l)i^«(*^)(/), (24) 

m— — 2 


which in turn implies 


ri2(/) — C'2ri2,2(/) 

2 

= C2 E ^™(7^*l,7^Ul)i?2t2™)(/)- 

m— — 2 


(25) 


Thus, the only non-zero contribution to the overlap re¬ 
duction function comes from the quadrupole gradient 
terms. Using Eq. (20), it follows that 


R2(2m)if) ~ Fm(JZU2,'TZV2) 


jo(a) + (-l)™+i(2-2 + l)2 

-h(-l)™(2-2-hl)(2-4-hl) 


f ^ 

2 

2 \ 

[2 

2 2 \ 

V —m 

m 

0 j 

1 2 

-2 0 j 


2 2 4 

—m m 0 


J2(c 


2 2 4 
2-2 0 


J4(a) 


(26) 


where a = 27r/|Aa;|/c. Thus, 


ri2(/) — C 2 [4ljo(a;) -I- Bj 2 {a) + Cji{a)\ , 


(27) 


where 


2 

E Fm{Tlui,TZvi)F^{TZu2,TZv2), 

m— — 2 
2 

5= E ^m(7^Ul,7^^l)E^(7^W2,7^^)2)(- 

m—— 2 

5 ^ ^ 

= - E F,n{'R-ui,TZvi)F^{TZu2,Tlv2){-^ 

m— — 2 
2 

C'= E i"m(7^Ul,7^^l)F;,(7^U2,771)2)(-!)" 

m—— 2 
12 . ^ ^ 

= — E F^{TZui,TZvi)F^{TZu2,TZv2)-7^ 

m =-2 ' 


Here we have used the definition of the Wigner 3-j sycalol 
to simplify the expressions for B and C. 

For an unpolarised, isotropic and uncorrelated back¬ 
ground, the above expression for ri 2 (/) in terms of spher¬ 
ical Bessel functions is similar in form to expressions for 
the overlap reduction function in Refs. [53] and (53] . The 
difference is that in those papers the expansion is in terms 


+ 2 2 2\ / 2 2 2 \ 

^ ^ \-m m 0 J \ 2 -2 0 J ’ 

2>^+1(to2-2), (28) 

m)!(2 -I- m)! ' 


of jo(a), ji{a)/a, and j 2 {a)fa'^, while the above expan¬ 
sion is in terms of joioi), j 2 (<a), and j 4 {a). These two 
expansions can be related using the recurrence relation 


2/ -I- 1 

jz+i(a) = ^jK«)-i/-i(«), (29) 

a 
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FIG. 1: Overlap reduction function ri 2 (/) for an unpolar¬ 
ized statistically isotropic background for the LIGO Hanford- 
LIGO Livingston detector pair, plotted on a logarithmic fre¬ 
quency axis. 


for which 


ri2(/) = C2 iA-B + C)jo{a) 

+ {3B - IOC) + 35C 


(30) 


A plot of the corresponding overlap reduction function 
for an unpolarised statistically isotropic background is 
shown in Fig. for the LIGO Hanford-LIGO Livingston 
detector pair. 


V. ROTATIONAL AND ORBITAL MOTION 


As mentioned in Sec. |Tj previous analyses for 
anisotropic stochastic backgrounds using ground-based 
interferometers have been formulated in terms of the 
distribution of gravitational-wave power on the sky li- 
[5]. In addition, these approaches typically use cross¬ 
correlation measurements between pairs of detectors as 
the input data. As such, these analyses are insensitive to 
the phase of the gravitational-wave background at differ¬ 
ent spatial locations. Recall that the position-dependent 
phase information appears explicitly in the Fourier com¬ 
ponents of the metric perturbations, 

(cf. Eq. ©)• It also appears in the response functions 
R^{f,k) and cf. Eqs. (15) and (20), where xq 

is the location of the detector at the time t at which 
the measurement is made. For cross-correlation mea¬ 
surements between a pair of detectors, the correlated re¬ 
sponse depends on the location of the detectors xi and 


X 2 only via their difference Ax = X 2 — xi, which is inde¬ 
pendent of choice of origin. Thus, the correlated response 
repeats after one sidereal day of observation. This means 
that the orbital motion of the Earth is effectively irrele¬ 
vant for such analyses; only the rotational motion counts. 
In fact, one can fold several days of observed data to a 
single sidereal day, and perform the analysis on the folded 
data [13 . This has obvious benefits in regards to the re¬ 
duction of data volume and the computational cost of the 
analysis. 

In contrast, the goal of our analysis is to recover both 
the amplitude and phase of the gravitational-wave back¬ 
ground at each point on the sky, based on a likelihood 
function that is not tied to cross-correlated data. As we 
shall see below, for such an analysis, the spatial loca¬ 
tions of the detectors are as important as their relative 
orientations. Since our detectors are ground-based inter¬ 
ferometers that orbit the Sun with the Earth, it is natural 
to reference our response functions and reconstructed sky 
maps back to the SSB. As such, we will take the origin 
of coordinates for our calculations to be the SSB. The 
detector locations are thus referenced from there. 

Due to the Earth’s rotational and orbital motion 
around the Sun, a single interferometer actually defines 
a set of virtual interferometers located along a quasi¬ 
circular ring 1 AU from the SSB. The Doppler shift in the 
observed frequency due to Earth’s velocity with respect 
to the SSB is not important for searches for broad-band 
gravitational-wave backgrounds, since x/c ^ 10“'^ intro¬ 
duces frequency shifts of at most Sf ^ few x 10“^ Hz 
in the frequency band relevant for ground-based interfer¬ 
ometers. Thus, the rotational and orbital motion of the 
Earth synthesizes a set of static virtual interferometers, 
each observing the gravitational-wave background from a 
different spatial location and with a different orientation. 

To compare the effects of rotational and orbital mo¬ 
tion, we calculate the time-scale over which measure¬ 
ments made by different virtual interferometers are cor¬ 
related. The relevant quantity is a = 27r/|Ax|/c, which 
appears in expressions for the overlap reduction func¬ 
tion Fi 2 (/), e.g. (27). But here Ax = xo(f 2 ) — ®o(D) 


is the spatial separation between the vertices of two vir¬ 
tual interferometers, defined by the vertex of a single 
(real) interferometer at times ti and t 2 ', and / is the 
frequency of a gravitational wave. The above relation 
can be turned into a correlation time-scale by writing 
|Ax| = vAt, where v is the average speed of the inter¬ 
ferometer (due to Earth’s rotational or orbital motion) 
over the time interval At = t 2 — ti, and then finding that 
value of At for which a = tt: 


2vf 


(31) 


This corresponds to a spatial separation |Ax| = A/2, 
where A = c// is the wavelength of the gravitational 
wave. For durations At < tcom measurements taken by 
the two virtual interferometers are correlated; for dura¬ 
tions At > tcorr, the measurements are uncorrelated with 
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one another. This is justified by noting that two detec¬ 
tors will (on-average) be driven in coincidence by a grav¬ 
itational wave propagating along their separation vector 
whenever its wavelength is more than twice the separa¬ 
tion between the detectors. This argument [51] provides 
a rough lower bound on the decorrelation timescale of 
the detectors, which will actually be slightly larger since 
we must average over all propagation directions of the 
gravitational waves when considering a stochastic back¬ 
ground. 

For a gravitational wave with frequency / = 100 Hz 
(A = 3 X 10® m) and v = 27ri?£;/(l day) « 500 m/s, 
which is relevant for Earth’s daily rotational motion, we 
find 


3000 s (rotational motion). 


(32) 


For V = 2ttRes/{^ yr) « 3 x 10^ rn/s, which is relevant 
for Earth’s yearly orbital motion, we find 


tn 


50 s (orbital motion). 


(33) 


Figure is the overlap reduction function ri 2 (/) for 
an unpolarized isotropic background, evaluated at / = 
100 Hz, for two virtual interferometers as a function of 
time. The left panel is for a set of virtual interferometers 
synthesized by the daily rotation of a detector positioned 
at the Earth’s equator, with no orbital motion. One can 
see that the detector decorrelates on a timescale of ~ 1 
hour, and recorrelates after 24 hours whenever it returns 
to its starting position. If we switch off daily rotation 
and synthesize a set of virtual interferometers from the 
orbital motion of the Earth around the Sun, then we get 
the overlap reduction function in the right panel. Since 
the orbital velocity of the Earth around the Sun is much 
larger than the velocity of a detector on the surface of the 
Earth, the virtual interferometers build up a larger sepa¬ 
ration baseline on a shorter timescale. Hence the overlap 
reduction function goes to zero much more rapidly in this 
case and will only recorrelate after 1 year. 

We investigate this decorrelation timescale by numeri¬ 
cally computing the overlap reduction functions for daily 
rotation and orbital motion at a variety of gravitational- 
wave frequencies. The times at which the detectors first 
decorrelate are shown in Fig. with the lower bounds 
given by Eq. (311. The decorrelation timescale does 


indeed obey a simple 1// scaling, and at 100 Hz the 
timescale for daily rotation and orbital motion are actu¬ 
ally ^ 67 min (= 4020 s) and 60 s, respectively. There¬ 
fore orbital motion of the Earth around the Sun will 
rapidly synthesize a large network of independent vir¬ 
tual interferometers from the motion of a single detector, 
with a resolving power that increases on a relatively short 
timescale. 


VI. MAP MAKING 

In this section, we extend the map-making formalism 
of [U [5] to data taken by a network of ground-based 


interferometers. The key observation is that the time- 
dependent ground-based interferometer analysis can be 
mapped to a static PTA-like analysis, for a set of static 
virtual interferometers in a quasi-circular ring 1 AU from 
the SSB. Unlike the static PTA analysis mm, the vir¬ 
tual interferometers are not all centered at the same lo¬ 
cation, but see the sky from different locations due to 
the Earth’s rotational and orbital motion. This allows 
for recovery of both the grad and curl components of the 
background, as discussed in Sec. |H|| in terms of the re¬ 
sponse functions. We shall demonstrate this explicitly 
via maximum-likelihood recovered sky maps in Sec. |VI E| 
below. 


A. Response vector 


As described in Sec. |HI[ the Fourier-domain response 
of detector / to a gravitational-wave background is 


g{/) = / or 

Js^ A 

(Im) P 


(34) 


where the response functions Rf{f,k), A 
P = {G,C} are given by Eqs. 
We write this response abstractly as 


= { -!-, x} and 
(jlbl) and (201. 


r = Rh, 


(35) 


where h denotes the components of the gravitational- 
wave background in either the pixel or spherical harmonic 
basis, and R denotes the corresponding response function 
in that basis. The response function R acts on h via a 
sum over polarizations and an integration over the sky, 
or a sum over polarizations and a sum over spherical 
harmonic components. 

When performing the data analysis to produce maps 
of the gravitational-wave sky, we need to discretize both 
the map (in terms of pixels or spherical harmonic com¬ 
ponents) and the observed data. This leads to a time- 
frequency decomposition where the data are broken up 
into segments of duration r, which should be short com¬ 
pared to the timescale over which the orientation of the 
detectors changes appreciably. Since the peak sensitivity 
of the advanced ground-based interferometers is ~ 100 
Hz, we take the minimum segment duration to be the 
time required for a detector to decorrelate from itself 
under orbital motion, thus synthesizing an independent 
virtual detector (e.g., r « 60 s). The longest segment du¬ 
ration is the time beyond which the Earth’s rotation will 
have appreciably changed the antenna response pattern 
orientation (e.g., r « 2048 s). Each segment of data is 
then discrete Fourier transformed, yielding a finite num¬ 
ber of components for the vectors r and h. In the follow¬ 
ing, we will denote the discrete (positive) frequencies by 
/j, where j = 1, 2, • • • , Nf] the sky pixels by /c„, where 
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FIG. 2: Overlap reduction function at / = 100 Hz for two virtual interferometers as a function of time. The left-hand plot 
is for a set of virtual interferometers located on Earth’s equator, associated with Earth’s daily rotational motion. The virtual 
interferometers have one arm pointing North and the other pointing East. There is no orbital motion for this case, as the center 
of the Earth is fixed at the SSB. The right-hand plot is for a set of virtual interferometers at 1 AU from the SSB, associated 
with Earth’s yearly orbital motion. There is no rotational motion for this case, as the interferometers are located at the center 
of the Earth in its orbit around the Sun, with the orientation of the interferometer arms unchanged by the orbital motion. 




FIG. 3: Analysis of Fig. [^repeated for 10 < / < 1000 Hz. The plots show the first time the overlap reduction function between 
virtual interferometers goes to zero for daily rotation {left panel) and orbital motion {right panel). The solid line in the left 
panel does not extend fully to 10 Hz since the overlap reduction function does not go to zero at low frequencies. A rough 
indication of when the overlap reduction function should go to zero is given by considering that a pair of detectors should 
be driven in coincidence by a passing gravitational wave when the wavelength is more than twice the separation between the 
detectors. This defines a lower bound on the decorrelation timescale of the virtual interferometers, shown by a dashed line. 


n = 1,2,--- ,fVpix; and the spherical harmonic compo¬ 
nents of the sky by {Im), where I = 0,1, • • • , Zmax, and 
The detectors (interferometers) are labelled 
by the index I = 1,2,-•• ,Nd, and the time segments 
recorded by detector I as tn, where i = 1,2, ••• ,Ni. 


Combining the responses from all detectors, we have 
h.={h+{fj,kn),h^{fj,kn)} or {a^^)(/,-)} , 

R ={R^^{fj, kn), Riiifj, kn)} 

or {Rniim)ifj),Rn(im)ifj)} ■ 

(36) 

Written this way, the time-dependent ground-based in¬ 
terferometer analysis is mapped to a static PTA-like anal¬ 
ysis, for a set of virtual interferometers synthesized by the 
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Earth’s rotational and orbital motion around the Sun. 
The response vector r has N = Nf complex com¬ 

ponents. In contrast to a cross-correlation analysis, we 
do not require that all the detectors have data for the 
same time periods. When the number of time segments 
Nt is the same for all Nj. detectors, then N = NtNdNf. 
Similarly, h is a complex-valued vector of dimension 
M = 2N-pi^Nf (pixel basis) or M = 2NmNf (spherical 
harmonic basis), where = (^max + 1)^ — 4 is the num¬ 
ber of spherical harmonic (Im) modes out to ^max- (The 
—4 in the last expression is because summations over I 
start at ^ = 2). The response function R is thus repre¬ 
sented by a complex-valued matrix of dimension N x M. 
Since the frequency components of r and h are identi¬ 
cal, the frequency transformation part of R is simply the 
identity matrix 

To simplify the discussion for the remainder of this 
section, we will work in the pixel basis. The subsequent 
calculations are formally identical in both bases, and the 
resulting sky maps are effectively the same, provided Zmax 
is chosen so the the total number of modes in the 
spherical harmonic basis is of the same order as the num¬ 
ber of pixels fVpix- 


B. Maximum-likelihood estimation 

Using the above notation, the measured data can be 
represented by an TV-dimensional complex vector d = 
{diiifj)}, with contributions, in general, from both the 
gravitational-wave signal r and detector noise n: 

d = r -I- n = Rh + n . (37) 

If we assume that the noise is Gaussian-stationary, then 
we can represent it by an iV x (Hermitian, positive def¬ 
inite) covariance matrix C, whose components are given 

by 


which is a multivariate Gaussian distribution for the 
noise. Note that there is no factor of 1/2 in the expo¬ 
nential as the matrix sum is over only positive-frequency 
components. Given the likelihood function, we can now 
use either Bayesian inference or frequentist maximum- 
likehood statistics to estimate the model parameters. 
The latter is relatively simple to do if we fix the noise, 
since the signal parameters enter linearly in the likelihood 
in Eq. (40). 


Maximizing the likelihood with respect to h leads to 


Hml = (R^C-iR)-iR^C-id 


(41) 


for the recovered map. This is only a formal expression, 
however, since the Fisher matrix. 


F = R^C-iR, 


(42) 


is not invertible in general, since R may have not have full 
column rank. This occurs if the number of data points N 
is less than the number of modes M that we are trying 
to recover, or if the response matrix has null (or nearly 
null) directions—i.e., particular gravitational-wave skies 
have hnuii to which the network of detectors is effectively 
blind. This is discussed further in Sec. VI G[ Thus, cal¬ 
culating hjviL will, in general, require some form of regu¬ 
larization m- 

To do the regularization, it is simplest to work with the 
whitened data d = Ltd and whitened response matrix 
R = LtR, where L is a lower triangular matrix defined 
by the Cholesky decomposition of the inverse covariance 
matrix, = LLt. (An alternative approach, based 

on the unwhitened response matrix R, is described in 
App. [B|) In terms of the whitened quantities, we have 
F = RTR and 


Lml = (RtR)-iRtd = R+d, 


(43) 


where 


(38) 

where Sf is the frequency resolution. (Past analyses for 
stochastic backgrounds using ground-based interferome¬ 
ters have typically used Sf = 0.25 Hz, which is much 
greater than the 1/r ~ 0.001 Hz frequency resolution 
associated with the duration of the short-term Fourier 
transform, see e.g., H [2S]. This amounts to working 
with a coarse-grained frequency series, obtained by av¬ 
eraging over neighboring frequency bins.) If we further 
assume that the noise is uncorrelated between different 
detectors, then 

Cir{f)=SirSi{f), (39) 

where Si{f) is the (one-sided) power spectral density of 
the noise in detector /. In terms of these quantities, the 
likelihood function for the data is 

p(d|C, h) cx exp [-(d - Rh)^C-i(d - Rh)] , (40) 


R+ = (R^R)-iR^ (44) 

is the so-called pseudo-inverse of R. As before, this is 
just a formal expression as the M x M matrix RtR is 
not invertible in general. However, it is always possible 
to define the pseudo-inverse Rt’ in terms of the singular 
value decomposition (SVD) of R: 

R = USV^, (45) 

where U and V are TV x TV and M x M unitary matri¬ 
ces, and S is an TV X TU rectangular matrix with (real, 
non-negative) singular values (ffc along the diagonal, and 
zeros everywhere else. (Without loss of generality, we 
can assume that the singular values are arranged from 
largest to smallest along the diagonal.) Then 

R+ = VS+U^, (46) 

where S''" is dehned by taking the reciprocal of each non¬ 
zero singular value of S, leaving the zeros in place, and 
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then transposing the matrix. In terms of the SVD of R, 
the maximum-likelihood estimator can be written as 

Hml = R+d = VS+U^d. (47) 

The expected value and variance of Hml are given by 
(hMu) = R^Rh, 

uar(hML) = (hMLhj^L) “ (hML)(h[iL) = R^(R’^)''', 

(48) 

where the expression for the variance assumes that the 
gravitational-wave signal is weak compared to the detec¬ 
tor noise. Although not explicit in the last expression, 
the variance of Hml does depend on the noise C, since 
R = LtRand C ^ = LL^. 

If the non-zero singular values of S vary over several 
orders of magnitude, it is usually necessary to set to zero 
(by hand) all non-zero singular values less than or equal 
to some minimum value dmin (e.g., 10“® times the largest 
non-zero singular value). This reduces the noise in the 
maximum-likelihood reconstruction, which is dominated 
by those modes that we are least sensitive to. So in what 
follows, when we speak of non-zero singular values of S, 
we will actually mean the singular values satisfying 

^ dniin- 


C. Sky map basis vectors 


Expression (471 for the maximum-likelihood estimate 


has several nice geometrical properties [2]. In particu¬ 
lar, the columns of U and V corresponding to the non¬ 
zero singular values of S have the interpretation of re¬ 
sponse range vectors and sky map basis vectors respec¬ 
tively, in terms of which the measured response Rh and 
the maximum-likelihood estimate Hml can be written as 
linear combinations. To see this, let U(j.) and V(j.) denote 
the fcth columns of U and V, and let r < min(A^, M) 
be the number of non-zero singular values of S. Then it 
follows from Eqs. (45) and (47) that 


Rh = ^crfe(v(fc) • h)u(fe) , 




hML = X! ^(u(fc) ■ d)v(fc), 


(49) 


k=l 


where dot product of two vectors a and b is defined by 
a • b = a^b. If we expand d = Rh -|- h, then we can also 
write 


hML = '^{'^(k) ■ h)v(fe) -t R+n. (50) 

k=l 

Note that this last expression for Hml involves the projec¬ 
tions of h onto for only the non-zero singular values 
of S. 


It is important to discuss in some detail those cases 
where there are fewer data points than modes we are try¬ 
ing to recover (i.e., N < M), or if there are certain modes 
of the gravitational-wave background that our response 
matrix is insensitive to. For either of these two cases, 
the system of equations d = Rh is under-determined, 
which implies that there exist multiple solutions for the 
recovered gravitational-wave background: 

h = R+d -h (Imxm - R^R)harb , (51) 


where harb represents an arbitrary gravitational-wave 
background. The particular solution that we have chosen 
for hML (given by Eq. (47) or (50)) ignores the term pro¬ 


portional to harb, setting to zero those modes that we are 
insensitive to. Our solution also sets to zero the variance 
of these modes, as can be seen from the expression for 
uar(hML) given in Eq. (48): 


uar(hML) = R+(R+)^ = VS+(S+)^V^ 


(52) 


which can be diagonalized by a similarity transformation 
involving Vb This yields S+(S+)f, which has M — r 
zeros along its diagonal. 

In a Bayesian formulation of the problem, things will 
be different, however, as one must also specify prior prob¬ 
ability distributions for the signal parameters, in addition 
to the likelihood function ( [40| ). For a signal parameter 
(or a combination of signal parameters) corresponding to 
a mode of the background that the detectors are insen¬ 
sitive to, the marginalized posterior will simply recover 
the prior distribution on this parameter (or combination 
of parameters), since the data are completely uniforma- 
tive about this mode. This is more in line with what we 
would expect for a mode that is unconstrained by the 
data, but such an analysis requires the specification of 
prior probability distributions which frequentist estima¬ 
tors, like Hml, do not provide. We therefore choose to 
construct our maximum-likelihood estimator such that it 
sets the modes that we are insensitive to equal to zero, 
and acknowledge the fact that we cannot say anything 
about them with our experiment. 


D. Sky maps, uncertainty maps, and SNRs 


To construct sky maps from the maximum-likelihood 
estimator hML, we need to either restrict attention to 
a particular discrete frequency fj or perform an aver¬ 
age over the different frequency components. In either 
case, the dimensionality of Hml reduces to 2iVpix, corre¬ 
sponding to the -f and x components of the estimated 
gravitational-wave background at each pixel on the sky. 
Uncertainty maps for these estimates are given by the 
square-root of the diagonal elements of the variance esti¬ 
mate given in (48), 


ctml = V diag [uar(hML)] • 


(53) 
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Similarly, we can construct signal-to-noise ratio (SNR) 
maps by simply dividing the estimates of and at 
each pixel on the sky by the corresponding values of ctml- 
Examples of such maps are given in Sec. |VI E[ 

E. Simulations 

We now illustrate the mapping procedure described 
above by constructing maximum-likelihood estimates of 
the real and imaginary parts of hj^{f,k) and h^{f,k) 
for three different simulated gravitational-wave back¬ 
grounds: a point source and a spatially-extended back¬ 
grounds having only gradient or curl modes. For simplic¬ 
ity, we consider only a single frequency component / = 
100 Hz, and we pixelize the sky using a HEALPix [35] 
grid containing A^pix = 768 pixels. (The sky map vectors 
Hml and hinj thus have dimension M = 2Npi^ = 1536.) 
We will work primarily with a network of Nd = 6 de¬ 
tectors, comprising both the existing and planned large- 
scale, ground-based laser interferometers LIGO-Hanford, 
LIGO-Livingston, Virgo, KAGRA, INDIGO, and AIGO. 
(Relevant information for each interferometer is given 
in Table |lj which is adapted from |5S|.) For compari¬ 
son, we will also consider a reduced network having just 
Nd = 3 detectors (LIGO-Hanford, LIGO-Livingston, and 
Virgo), which is more realistic for the near future. The 
measured data will be given in the frequency domain, 
corresponding to short-term Fourier transforms of time 
segments of duration t = (1 sidereal day)/60 « 1436 s. 
The simulations for the 6-detector network will have a 
total of Nt = 400 samples for each interferometer, corre¬ 
sponding to 6.67 days of simulated data. The simulations 
for the 3-detector network will have either Nt = 400 or 
Nt = 800 samples for each interferometer, corresponding 
to either 6.67 days or 13.33 days of simulated data. The 
data and response vectors d and r will thus have dimen¬ 
sions N = NdNt = 2400 for the 6-detector network, and 
N = 1200 or 2400 for the two 3-detector networks. 

The detector noise will be described by an V x V 
block-diagonal covariance matrix, whose Nd blocks (cor¬ 
responding to the Nd detectors in the network) are each 
proportional to the unit matrix IjVjxAff The propor¬ 
tionality constants are the values of the one-sided power 
spectral densities £'/(/), I = 1,2, - ■■ ,Nd evaluated at 
/ = 100 Hz (see the last column of Table |T]) divided by 
4(5/, where Sf is the size of the frequency bins. The 
factor of 4 is due to the use of one-sided power spec¬ 
tral densities (one factor of 2) and the summation over 
only positive-frequency bins (the other factor of 2). For 
our simulations, we take Sf = 0.25 Hz, as is common 
for stochastic background searches using ground-based 
interferometers [26| . The real and imaginary parts of 
the noise vector n are generated by randomly drawing 
independent samples from a multivariate Gaussian dis¬ 
tribution defined by this block-diagonal matrix. 

The simulated gravitational-wave backgrounds will 
consist of a point source and two spatially-extended 


distributions. The point source is not an ideal point 
source, but more of a Gaussian ‘blob’, since we gen¬ 
erate it out to only = 10 (see the first row of 

maps in Fig. [^. Nonetheless, it serves its purpose as 
being a simple yet extreme example of an anisotropic 
background for both h+ and hx- (We also considered 
single-pixel point sources and found similar results, but 
the corresponding sky maps are not very clear.) The 
two spatially-extended backgrounds are a grad-only sta¬ 
tistically isotropic background with equal contributions 
for multipoles 2 < ^ < 10, and a curl-only statistically 
isotropic background with equal contributions for multi¬ 
poles 2 < I < 5 (see the second and third rows of Fig. ® . 
These last two backgrounds were also considered in ^ 
in the context of pulsar timing arrays. It is interest¬ 
ing to compare the recovered sky maps for the ground- 
based and pulsar timing analyses, especially for the curl- 
only background, which cannot be recovered using timing 
residual data from a pulsar timing array. The amplitudes 
of the injected gravitational-wave backgrounds were cho¬ 
sen to give reasonable recoveries after just a few days of 
simulated data. For the values of the noise spectral den¬ 
sities Si{f) given in Table [ij we found that an amplitude 
A = 4x 10“^® was sufficient for the three different back¬ 
grounds. (If we used a smaller value of A, we would have 
had to integrate for a longer period of time.) The faithful¬ 
ness of the recovery is measured by calculating the match 
between the injected and maximum-likelihood-recovered 
sky maps, 

2 (^inj ' h]V[L “f hiy[L ’ hjnj) 

M = — , - - . (54) 

Y hinj ' hinj V hML ' h]V[L 

This is just the coherence between the two maps. We 
will also construct uncertainty maps and SNR maps to 
evaluate how well we can recover the injections. 

Figure[^is a plot of the match as a function of the num¬ 
ber of days of observation for the 6-detector network and 
a noiseless point-source injection. The match increases 
as the total observation time increases as expected. Note 
we get perfect match after 5 days of observation. This 
follows from the fact that the total number of data points 
taken by the 6-detector network over 5 days is given by 
N = 5 days X 60 samples/day x 6 = 1800, which is greater 
than the number of modes M = 2Npix = 1536 we are 
trying to recover. Thus, in the absence of noise we have 
(more than) enough information to completely recover 
the injected background after 5 days of observation. We 
would have complete recovery for the two other simulated 
backgrounds as well. 

Sky maps of the recovered point-source background in¬ 
jected into noisy data are shown in Fig. [^ These maps 
are for the 6-detector network with N = 2400 total data 
points, corresponding to 6.67 days of total observation. 
The first row shows the injected background. The second 
row shows the maximum-likelihood sky map estimates, 
which are the real and imaginary parts of the /i+, hx 
components of Hml- The third row shows the uncertainty 
maps, as specified by ctml, and the fourth row shows the 







12 


Interferometer 

Longitude 


Latitude 

Orientation 

Spectral density 

LIGO, Hanford 

119° 24' 27.6" 

’ W 46° 27' 18.5" N 

279.0° 

1.591 X 10"'‘' Hz-" 

LIGO, Livingston 

90° 46' 27.3" 

W 

30° 33' 46.4" N 

208.0° 

1.591 X 10"'^'^ Hz"^ 

Virgo, Italy 

10° 30' 16" 

E 

43° 37' 53" N 

333.5° 

2.063 X 10"'^'^ Hz-1 

KAGRA, Japan 

137° 10' 48" 

E 

36° 15' 00" N 

20.0° 

9.320 X 10-1® Hz-1 

INDIGO, India 

74° 02' 59" 

E 

19° 05' 47 N 

270.0° 

1.591 X 10-11 Hz-1 

AIGO, Australia 

115° 42' 51" 

E 

31° 21' 29" S 

45.0° 

1.591 X 10-11 Hz-1 


TABLE I: Geographic information for ground-based interferometers used in our simulations, adapted from [28] . Orientation is 
the angle that the bisector of the two interferometer arms makes with geographic North (positive for directions pointing East 
of North). All interferometers are assumed to have 90° opening angle between the two arms. Spectral density is the value of 
the one-sided noise power spectrum S„(/) for the corresponding interferometer, evaluated at / = 100 Hz. The values of 5„(/) 
for LIGO, Virgo, and KAGRA are taken from design sensitivity documents and publicly accessible data | 29H31| : the values for 
INDIGO and AIGO are taken to be the same as those for the LIGO interferometers, as they are in the initial planning stages. 
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FIG. 4: Simulated maps (at a single fixed frequency) for three different anisotropic gravitational-wave backgrounds: (i) point 
source located at 40° latitude, 60° longitude having /max = 10 (first row); (ii) grad-only statistically-isotropic background with 
Cl = const for 2 < / < 10 (second row); (iii) curl-only statistically-isotropic background with Ci = const for 2 < Z < 5. The 
four columns correspond to the real and imaginary parts of /i+ and hx- 


SNR maps. The max SNR at the location of the point 
source is approximately 10. The match is /i = 0.64 for 
this particular simulation. 

Note that the uncertainty maps for the real and imag¬ 
inary parts of (or h^) are the same. The uncertainty 
values are also fairly constant over the sky, with val¬ 
ues around 3 x 10“^®. Thus, the SNR maps look very 
similar to the maximum-likehood maps but, of course, 
have different values since they represent different quan¬ 
tities. It is also the case that the uncertainty maps for 
the other simulated backgrounds (grad-only and curl- 
only background) will be identical to that for the point 
source background, since Eq. ( [5^ for (Tml depends only 
on the response matrix R and noise covariance matrix C 


via (481—i.e., it is independent of the background that 
one is trying to recover, at least in the weak-signal limit, 
which we have assumed in our analyses. So for the other 
two simulated backgrounds, we will show only the in¬ 
jected and maximum-likelihood recovered sky maps, and 
not the uncertainty and SNR maps. 


Figure [7] is identical to Fig. but for the 3-detector 
network having the same number of total data points 
{N = 2400) as the 6-detector network. The total obser¬ 
vation time is thus twice as long, in order to compensate 
for the reduction in the number of interferometers. Note 
that the uncertainty maps have values that are slightly 
larger than for the 6-detector network. Also, the match 
is p- = 0.59, which is slightly smaller than that for the 
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FIG. 5: Match as a function of number of days of observa¬ 
tion for the 6-detector network and a noiseless point-source 
injection. A match value equal to 1 corresponds to perfect 
recovery. 


6-detector network. Thus, we see from this simulation 
that the 3-detector network performs nearly as well as 
the 6-detector network if we integrate long enough to 
acquire the same number of total data points. The per¬ 
formance of the 3-detector network is much worse than 
the 6-detector network if we integrate for the same obser¬ 
vation time, since then the total number of data points 
for the 3-detector network is half as large (see the fourth 
rows of Figs. and below). 

Maximum-likelihood recovered maps for the grad-only 
and curl-only backgrounds using the 6-detector and two 
3-detector networks mentioned above (one having the 
same number of total data points as the 6-detector net¬ 
work; the other having half as many data points) are 
shown in Figs. and The corresponding match val¬ 
ues for the grad and curl recoveries are /i = 0.81 and 
0.85 for the 6-detector network. For the two 3-detector 
networks, the match values are fi = 0.80 and 0.84 when 
the total number of data points is the same as for the 6- 
detector network, and /r = 0.55 and 0.60 when the total 
number of data points is half as many. The SNR values 
as a function of sky location for the different recoveries 
range from about —5 to 5 for the strong recoveries and 
—4 to 4 for the weaker recoveries. As can be seen from 
the fourth row of these two figures, when the 3-detector 
network has only half as many total data points as the 
other two networks, the structure in the grad-only and 
curl-only sky maps is not nearly as clearly recovered as 
for the other detector networks. 

The most important take-home message is that the 
grad-only and curl-only backgrounds can both be recov¬ 
ered with a network of ground-based interferometers. 
This is in contrast to the case for a pulsar timing array. 


which is completely insensitive to a curl-only background 
(see [I], and in particular Fig. 11 from that paper). As 
mentioned earlier, the rotational and orbital motion of 
the Earth synthesizes a set of virtual interferometers that 
sample the gravitational-wave field from many different 
spatial locations. This allows for the reconstruction of 
both grad and curl modes of the background, unlike the 
case for pulsar timing arrays. 


F. Minimum duration between data segments for 
independent measurements 


Having shown our map recovery techniques to be suc¬ 
cessful in the context of noisy simulated injections, we 
now return to an issue discussed at the end of Sec. El- 
namely, the minimum time increment between observa¬ 
tions required to synthesize a network of independent vir¬ 
tual interferometers, and thus avoid degeneracies in the 
information content of our measured strain signal. We 
consider two cases: (i) a single AdvLIGO Livingston de¬ 
tector, and (ii) the full 6-detector advanced network pre¬ 
viously discussed. In both cases we assume a total of 
1200 strain measurements of the gravitational-wave sky 
have been recorded. As previously, the h+, compo¬ 
nents of the sky are decomposed into 768 pixels for a 
total of 1536 unknown parameters to be determined by 
our search. We compute maximum likelihood maps from 
the 1200 observations, which are carried out over various 
total timespans to investigate how the match of the re¬ 
covered map with the injected map scales with At, the 
time between observations. 


Our results are summarized in Fig. For all cases, 
we find that the match of the recovered maps with the in¬ 
jected map is poor for small time increments between ob¬ 
servations, since the detector(s) will not have moved far 
enough to establish independence from its previous posi¬ 
tion. With only orbital motion of a single detector, the 
match values are only able to plateau at ~0.5. Adding in 
the influence of daily Earth rotation seems to ameliorate 
this poor match behaviour. The origin of this effect can 
be deduced from the singular values of the response ma¬ 
trix in both cases. This is shown in Fig. 11 for a single 
LIGO Livingston detector, where the addition of Earth 
rotation acts to break degeneracies in the response matrix 
and conditions it to have a much smaller dynamic range 
of singular values. Rotating the Earth acts to sweep the 
antenna beam pattern of the detector across the sky and 
provides additional information with which to measure 
the gravitational-wave background. With only orbital 
motion the arms of the detector remain in fixed orien¬ 
tations, and hence so does the detector’s antenna beam 
pattern. We also show in Fig.[^the match behaviour for 
a full 6-detector advanced network. In this case we al¬ 
ready have information from multiple orientations of the 
antenna beam patterns by virtue of the different global 
placements of the detectors. Hence, the inclusion of daily 
Earth rotation makes little impact on the match value. 
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FIG. 6: Recovery of the simulated point source in noise for the 6-detector network. Injected maps (first row); maximum- 
likehood recovered maps (second row); uncertainty map (third row); SNR map (last row). Note that the uncertainties maps 
for the real and imaginary parts of (or hx) are the same. 


which plateaus at ~0.9 even with orbital-only motion. 

The key lessons here are that Earth’s daily rotation is 
an important influence on top of the orbital motion of the 
Earth around the Sun, since it sweeps the detector an¬ 
tenna beam patterns across the sky to gather additional 
information about any gravitational-wave signal of inter¬ 
est. Furthermore, from Fig.[T^we can clearly see that the 
first peak of the match value occurs at ^ 50 — 60 s, when 
the detectors decorrelate from themselves for the first 
time (see Fig. and are no longer driven in coincidence 
by a passing gravitational wave. With this time incre¬ 
ment the detector’s strain measurements are effectively 
independent from their preceding or subsequent measure¬ 
ments, thereby allowing us to synthesize a large network 
of virtual interferometers from the daily and orbital mo¬ 
tion of the Earth. The small dip after the first peak may 
be due to the detectors being driven in anti-coincidence, 
thereby losing some of their independence. However, the 
match value recovers in the limit of large At, since the de¬ 
tectors are then separated by several gravitational-wave 
wavelengths and this behaviour is averaged out. 


VII. DISCUSSION 


We have presented a new method for mapping the 
gravitational-wave sky using a network of ground-based 
laser interferometers. This method extends the for¬ 
malisms developed in mm. which were originally applied 
to the case of pulsar timing arrays. We have shown that 
we can recover both the gradient and curl components 
of a gravitational-wave background, as a consequence of 
the spatial separation of the individual interferometers 
in the network, or of a single interferometer at different 
times during its rotational and orbital motion around 
the Sun. This is in contrast to the case for a pulsar 
timing array, which is completely insensitive to the curl 
modes. Also, by mapping both the amplitude and phase 
of h^{f,k) and hx{f,k) as functions of direction on the 
sky (as referenced from the SSB), our method extends 
previous approaches [3H9] for anisotropic backgrounds, 
which map the distribution of gravitational-wave power, 
|/i+P -|- l^x P- Our formalism can be cast in terms of ei¬ 
ther the traditional -I- and x polarization modes of the 
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FIG. 7: Same as Fig. but for the 3-detector network having the same number of total data points (A^ = 2400) as the 
6-detector network. Injected maps (first row); maximum-likehood recovered maps (second row); uncertainty map (third row); 
SNR map (last row). The maps are more-or-less the same as for the 6-detector network shown in Fig. 


background {/i+(/, k), (/, fc)}, or the gradient and curl 

modes {a^^)(/),a^^^(/)}, with respect to a decomposi¬ 
tion of the metric perturbations in terms of spin-weighted 
or tensor (gradient and curl) spherical harmonics [T]. 

The results of the simulations presented in Sec. |VIE 
can be thought as a proof-of-principle demonstration of 
the general map-making formalism described in the rest 
of the paper. The actual analysis of real data from a net¬ 
work of advanced interferometers will most likely differ 
from this simplified scenario in several ways: 

(i) The amplitudes of the simulated backgrounds were 

chosen to be sufhcently large, so as to allow for fairly de¬ 
cent recovery after only a few days of observation. Much 
weaker backgrounds will require an increased observation 
time, of order months or years, noting that the (power) 
signal-to-noise ratio scales like where A is the am¬ 

plitude of the background and T is the total observation 
time. 

(ii) For initial analyses, it might be easier to work in 
the tensor spherical harmonic basis, and estimate the 
grad and curl components of the background {apm)(/), 


some relatively small value of Zmax, e.g., 
l-max = 10. This would reduce the number of modes that 
we would need to recover from 2A^pix (= 1536 for exam¬ 
ple) to 2iVniodes = 234 at each discrete frequency. The 
estimates of the grad and curl components can th en be 
converted to sky maps of h+{f,k), h^{f,k) using (12). 


(iii) Varying noise levels in the detectors (on a time 
scale > the segment duration r of our short-term Fourier 
transforms) will complicate somewhat the expression for 
the noise covaraince matrix C. The Nd block ma¬ 
trices that enter the expression for C will no longer 
be proportional to the unit matrix IjVjxWtJ tint rather 
will have diagonal elements proportional to Si{f;tji), 
i = 1, 2, • • • , Nt, reflecting the time-varying noise levels 
in detector I. 


(iv) As ground-based interferometers are broad-band 
detectors, we will have measurements at a set of discrete 
frequencies /j-, j = 1, 2, • • • ,Nf, where V/ ~ several hun¬ 
dred to a few thousand depending on the frequency bin 
size Sf. For initial analyses, it will probably be simplest 
to average the estimates of h+{fj,k), hx{fj,k) over the 
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FIG. 8: Recovery of the grad-only background in noise. Injected maps (first row); recovered maps for the 6-detector network 
(second row); recovered maps for the 3-detector network having the same number of total data points (N = 2400) as the 
6-detector network (third row); recovered maps for the 3-detector network having half as many total data points (N = 1200) 
as the 6-detector network (fourth row). 


different frequency components. 

(v) If one would like to compare the consistency of dif¬ 
ferent models of a stochastic background with the mea¬ 
sured data—e.g., is the measured data consistent with an 
unpolarized, isotropic background or with a background 
having a non-zero dipole component or with correlated 
emission on the sky, etc.—a Bayesian formulation of the 
problem would be more appropriate. The different mod¬ 
els would be defined by the appropriate choice of vari¬ 
ables for the stochastic background and prior probably 
distributions for these variables. Bayesian model selec¬ 
tion would then be used to select between the competing 
models. 

Perhaps the most compelling reason for using the for¬ 
malism presented here is that it provides a completely 
generic approach to mapping the gravitational-wave sky. 
It allows us to construct a map of the background that 
extracts all of the information that is possible to extract 
from the measured data. With the advanced ground- 
based interferometers coming on-line at the end of this 
year, and with the first detection of gravitational waves 


expected to follow shortly thereafter, it seems appropri¬ 
ate to utilize approaches such as this that attempt to 
maximize the science return of the data. 
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FIG. 9: Recovery of the curl-only background in noise. Injected maps (first row); recovered maps for the 6-detector network 
(second row); recovered maps for the 3-detector network having the same number of total data points {N = 2400) as the 
6-detector network (third row); recovered maps for the 3-detector network having half as many total data points {N = 1200) 
as the 6-detector network (fourth row). 
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Appendix A: Derivation of gradient and curl 
response functions 

Here we derive the gradient and curl response func¬ 
tions for an interferometer in the small-antenna limit, 


allowing for a non-zero displacement xq of the vertex of 
the interferometer from the origin of coordinates. 

Expressions for the response functions evaluated in a 
reference frame whose origin is located at the vertex of 
the interferometer were derived in Appendix D of 

Rflm)if) =Sl2Y\l\ [Y 2 m{u) - Y 2 miv)] , 

where m, v are unit vectors in the directions of the two 
arms of the interferometer. We have put bars on the 
above expressions to distinguish them from similar un¬ 
barred quantities that we will calculate in a reference 
frame whose origin is at the solar system barycentre 
(SSB). Note that is independent of frequency and 

is non-zero only for the quadrupole modes, I = 2. 

Under a translation of reference frames from the SSB 
to the vertex of the interferometer located at xq, the 
Fourier components habif-, k) of the metric perturbations 
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FIG. 10: Match values for various noisy map injections (averaged over 500 noise realizations) as a function of the time increment 
between 1200 observations, for a single detector and a 6-detector network. For a single detector, orbital-only motion is only able 
to achieve a match value of ~ 0.5 at large At. Daily rotation provides additional information by sweeping the antenna beam 
pattern across the sky, thus giving excellent plateau match values of ~ 0.9. The global placement of a network of detectors 
(and their differing antenna beam pattern orientations) gives excellent match values even for orbital-only motion, and daily 
rotation does not significantly improve this match. We see that the first peak in the match value occurs at ~ 50 — 60 seconds 
(shown as a grey strip), when the detectors first decorrelate from themselves and are no longer driven in coincidence by a 
passing gravitational wave. The small dip after the first peak is due to the virtual detectors being driven in anti-coincidence 
by the gravitational wave, thereby losing some of their independence and diminishing the match. However at larger At this 
behavior is averaged out over several gravitational-wave wavelengths, allowing the match value to recover. 
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FIG. 11: Singular values of the response matrix for a single LIGO-Livingston interferometer, for 1200 total observations with 
At = (1 sidereal day)/60 « 1436 s. (The singular values are normalized by the largest singular value, corresponding to the first 
diagonal element.) The dashed blue line shows the singular values when the detector is affixed to an Earth undergoing orbital 
motion only, whilst the solid blue line shows the singular values when the Earth is both rotating and orbiting. This indicates 
that in the orbiting-only case, regularization of the response matrix at machine-level precision (e ~ 10“^®) will not remove 
the very small singular values after diagonal element 768, requiring a more stringent cutoff level. In contrast, when the extra 
influence of the Earth’s rotation is introduced, the additional information provided by the changing detector arm orientations 
(which sweep the antenna beam pattern across the sky) acts to drastically improve the conditioning of the response matrix. 
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hab{t,x) in the “cosmic” (or SSB) frame transform to 

habU, k) = habif, fc)e-*2-/^--o/c (A2) 

in the detector frame. The correponding mode expan¬ 
sions in the two frames are given by 

hab{f,k)=Y.T.^(lm)if)YiL)abik), 

(Im) P 


Pb{f,k) = Y.Y.^ilrn)if)Y(L)abCk)- 

-- 


(A3) 


This last equation for hab{f, k) ca n be inverted to find 
®^m)(/) terms of using (A2) and the orthog¬ 

onality of the gradient and curl spherical harmonics: 




/S2 




/S2 


(A4) 


{I'm') P' 


IS^ 


Using the identity: 


g * 27 r/fc.:ro/c ^ ^ (io)yLM (fc) , O = 27r/|xo |/c , 


(AS) 


we obtain 


L=0 


oo L 


M^-L 


n 

dm)(/)= E E4'-')(/)E E M-^fjL{a)Yl^{xo) d2%y({);„,),,(ft)Y(L)“'>*(fc)yiM(fc) (A6) 

{I'm’) P' L^r\M=-r, 


L=0 M=-L 


relating the mode coefficients in the two frames. 

To make the connection between the mode coefficients and the corresponding response functions, we note that the 
detector response f{f) (or r{t)) to the gravitational-wave background will have the same value regardless of which 
frame we choose to evaluate it in. Thus, 


V) = J2J2Klm)4m)if) 


{Im) P 
2 




m— — 2 
2 


oo L 


= E ^f2’n) E E4'™')(/)E E 47r(-*)^J^(a)r;^(io) / y(E')a.(^)d^(L)“'’*(fc)d^LM(fc) 

m=-2 {I'm') P' L=0M=-L 

= E E<-')(/)4'-')(/)’ 

{I'm') P' 

(A7) 


where 


2 OO L 


R{lm){f)= E E E Rf 2 m')M-^fjLia)Ylj^ixo)[ d^^ U(L)a. (fc) ■ (AS) 

m' = -2 L=0 M=-L 


We can write this last expression explicitly in terms of Wigner-3j symbols if we replace the grad and curl spherical 
harmonics in the integral by spin-2 spherical harmonics using: 


^(L)ah(fc)^(E)“'*(fc) = i -2Ylm{k)-2Y^^,{k) + ^Y^ikhY*^, (k) 


1 r 


Y{km)ab(k)Y^2m’)^^*{k) = ^ -2Ylm{k)-2Y,*^, {k) - 2YUk)2Y*^, {k) 


rG ab*/ 


2 L 
1 r 


2f L 


(A9) 
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This leads to 

2 oo L 

Rflm)U)= H H H 

m'——2 L—0 M— — L 

(-l)’"' /(2-2 + l)(2Z + l)(2L+l) f 2 


Att 


I L 
n M 


m'——2 L—0 M— — L 

(-1)™' I{2 ■ 2 + 1){21 + 1){2L + 1) ( 2 


2i 


47r 


I L 
TL M 


2 I L 
-2 2 0 


2 I L 
-2 2 0 


+ 


2 I L 
2-2 0 


2 I L 
2-2 0 


(AlO) 


These expressions can be further simplified using the symmetry property 

_ ^_^^1+^2+^ 


h h I 

mi m 2 m 


h h I 

—mi —m 2 —m 


to eliminate one the Wigner 3-j symbols in terms of the other, and the triangle inequality 

1^1 — ^21 ^ ^ ^ “ i “ ^2 ^ — 2 , L I ‘2 

to collapse the infinite sums over L to sums over just 5 terms. The final expressions are 

2 /+2 L 

R{im)if)= H H H ^(2m')47r(-*)qL(a)nM(io) 

m'——2 L—l — 2 M——L 

(-1)"*' I{2 ■ 2 + 1){21 + 1){2L + 1) ( 2 


(All) 


(A12) 


47r 


I L 
n M 


2 I L 
2-2 0 


[(_l)'+i + 1] , 


2 Z+2 L 


(A13) 


m'——2 L=l — 2 M— — L 

(-l)’"' /(2-2 + l)(2Z + l)(2L+l) 


2i 


The Wigner 3-j symbol selection rule 


47r 


(2 I L \ (2 I L\ . . 

-to' m M ) \ 2 -2 Q ) > J ■ 


— to'+to + M = 0 


(A14) 


implies that the sum over M collapses to only those values satisfying M = m! — m and \M\ < L. Note also that in a 
reference frame with the z-axis chosen along xq, 




Y*j^(0,^) = Smo 


2L + 1 


4tt 


(A15) 


so for this case the sum over M reduces to just the M = 0 value. The selection rule —to' + to = 0 and \m'\ < 2 then 
imply non-zero values for only \m\ < 2 in this special frame. 


Appendix B: Equivalence of whitened and 
non-whitened analyses 

As discussed in Sec. IWBl we are interested in find¬ 
ing the maximum-likelihood value Hml of the likelihood 
function 

p(d|C,h)cxexp[-(d-Rh)lC-i(d-Rh)] , (Bl) 


when the Fisher matrix F = RlC“^R is not invertible. 
One approach, described in [T] , is to work with the SVD 
of the response matrix R: 

R = USVl. (B2) 

In general we can write U = [Ur.U„], where Ur- is an 
N X r matrix denoting the range of the response matrix 
R, where r equals the number of non-zero singular values 
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in S. We then replace Rh in the likelihood by U^b, 
where b is a vector of dimension r, and then proceed 
as for a non-singular response. The maximum-likelihood 
value for b is then 


and obtaining Hml from bML using Eq. (B4). The first 
approach requires no modification, but for the second ap¬ 
proach we now need the SVD of the whitened U^. matrix, 

LtU.: 


bML = (ui;c-iu,)-iu^c-M , (B3) 


= USV^, (B8) 


and the corresponding maximum-likelihood estimate of 
the gravitational-wave sky is 

bML = VS+Bml , (B4) 


where is the r x M dimensional matrix obtained by 
crossing out the last N — r rows of S, and is the 
pseudo-inverse of S,.) obtained by taking the reciprocal of 
each non-zero singular value of and then transposing 
the resulting matrix. 

A n alt ernative approach, which we described in 
Sec. VIB is to work with the the whitened data d = Lid 


and whitened response matrix R = L^R, where L is a 
lower triangular matrix defined by the Cholesky decom¬ 
position of the inverse covariance matrix, = LLb 
Working with the SVD of R: 


R = USVl 


(B5) 


we have 


Hmt = VS+Uld. 


(B6) 


Now the SVD of a product cannot be simply written in 
terms of the SVDs of the individual matrices. However, 
we can show the equivalence of these two approaches for 
maximizing the likelihood, by working with the equiva¬ 
lent likelihood that was introduced in the first approach, 
i.e., 

p(d|C,h) oc exp [-(d-U^b)lC-i(d-Urb)] , (B7) 


for which 


bML = vs+trid. 


(B9) 


Since we are now working only with the range of R and 
the noise covariance matrix C is positive definite, the 
rank of LlU^ must equal the rank of U,,. As before, we 
can write U = [Uf.U„], where U,. is an V x r matrix, 
which gives the range of LlUr. Thus, we can equivalently 
write Eq. (B8) as 


LlU,. = , (BIO) 


where Hr is an invertible, square r x r matrix obtained, 
as before, by crossing out the last N — r rows of S. Erom 
this last equation we now see that 

U,. = LlU^VS-^, (Bll) 

and 

bML = VS^^Uld 

= VS-^S-iVlUlLd (B12) 

= (uic~iu^)-iuic-M, 


where the final equality follows from the observation that 
UtC-^U^ = UlLLtU,. = VS^Vt (which is a conse¬ 
quence of Eq. ( BIO )). We have thus recovered the result 
given in Eq. (B3), which was obtained without whitening 
the data. 
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